##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
suppressWarnings(library(edgeR, quietly = T))
library(RColorBrewer)
library(ggplot2)
library(ggrepel)
library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)
library(ggsci)
library(ggrastr)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character"),
    make_option(c("--foldchange_t"), type = "character"),
    make_option(c("--q_t"), type = "character"),
    make_option(c("--subtype"), type = "character")
)

if(1!=1){
    

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file
q_t <- as.numeric(opt$q_t)
foldchange_t <- as.numeric(opt$foldchange_t)
subtype <- opt$subtype

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

class_type <- c( "IGC" , "DGC")
##########################################################################################

count_norm <- dat_tpm
rownames(count_norm) <- count_norm$gene_id
count_norm <- count_norm[,-which(colnames(count_norm) == "gene_id")]
##只提取全部或某个亚型的IGC和DGC样本
if(subtype=="All"){
  count_norm <- count_norm[,grep( "_IGC|_DGC", colnames(count_norm) , value = T)]
} else {
  use_sample <- info$Tumor[match(sapply(strsplit(colnames(count_norm), "_"), "[", 1), info$Tumor) & info$Molecular.subtype==subtype]
  count_norm <- count_norm[,grep( "_IGC|_DGC", colnames(count_norm) , value = T)]
  count_norm <- count_norm[,which(sapply(strsplit(colnames(count_norm),"_"),"[", 1) %in% use_sample)]
}


conditions <- factor(sapply(strsplit(colnames(count_norm) , "_") , "[" , 2))
table(conditions)
print(colnames(count_norm))
##########################################################################################
## 多组，两两做差异表达
result_diff <- c()

    class1="DGC"
    class2="IGC"

    ## Run the Wilcoxon rank-sum test for each gene
    count_norm_use <- count_norm

    pvalues <- sapply(1:nrow(count_norm_use),function(i){
        data<-cbind.data.frame(gene=as.numeric(t(count_norm_use[i,])),conditions)
        p <- wilcox.test(gene~conditions, data)$p.value
        return(p)
        })
    fdr=p.adjust(pvalues,method = "fdr")

    ## Calculate the fold-change for each gene
    conditionsLevel <- c(class1 , class2)
    dataCon1 <- data.frame(count_norm_use[,c(which(conditions==conditionsLevel[1]))])
    dataCon2 <- data.frame(count_norm_use[,c(which(conditions==conditionsLevel[2]))])
        if(ncol(dataCon1) > 1) {
            row_means1 <- rowMeans(dataCon1)
            } else if(ncol(dataCon1) == 1) {
                row_means1 <- dataCon1[[1]]
            }
        if(ncol(dataCon2) > 1) {
            row_means2 <- rowMeans(dataCon2)
            } else if(ncol(dataCon2) == 1) {
                row_means2 <- dataCon2[[1]]
            }
        foldChanges <- log2(row_means2/row_means1)
    ## Output results based on FDR threshold
    outRst<-data.frame(log2foldChange=foldChanges, pValues=pvalues, FDR=fdr)
    rownames(outRst)=rownames(count_norm_use)
    

        ## class1代表分子
    outRst$class1 <- class2
    outRst$class2 <- class1
    outRst$meanTMM_class1 <- row_means2
    outRst$meanTMM_class2 <- row_means1
    outRst$gene_id <- rownames(outRst)
    outRst=na.omit(outRst)
    
    result_diff <- rbind( result_diff , outRst)


##########################################################################################

colnames(result_diff) <- c("log2FoldChange" , "pvalue" , "padj" , "class1" , "class2" , "meanTMM_class1" , "meanTMM_class2" , "gene_id")

image_name <- paste0( out_path , "/DiffGene_IGC_DGC_",subtype,".tsv" )
write.table( result_diff , image_name , row.names = F , sep = "\t" , quote = F )

##########################################################################################
dat_diff <- result_diff
dat_diff <- merge( dat_diff , dat_gtf , by = "gene_id" )
##########################################################################################
diffplot <- function(use_dat = use_dat , image_name = image_name,subtype=subtype,Class1,Class2){

    #对原数据进行处理
    use_dat$log10_q_Value <- -log10(use_dat$padj)
    use_dat$log_FC <- use_dat$log2FoldChange
    use_dat$gene <- use_dat$Hugo_Symbol

    data <- use_dat[,c('gene','log_FC','log10_q_Value')]
    colnames(data) <- c('gene','log_FC','-log10_P_Value')

    #设置阈值
    logFC_cutoff <- log2(foldchange_t)
    log10_P_Value_cutoff <- -log10(q_t)

    options(ggrepel.max.overlaps = 20)

    plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
      geom_point_rast(data = subset(data,abs(log_FC)<logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point_rast(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point_rast(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'red',alpha = 0.4)+
      geom_point_rast(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'darkgreen',alpha = 0.4)+
      theme_bw()+
      theme(legend.title = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = 'none',
            axis.line = element_line(colour = "black"))+
      labs(x='log2(fold change)',y='-log10(adjusted p-value)')+
      geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
      geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
      geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                            aes(label = gene),size = 3,col = 'black')+
      ggtitle(paste0(subtype," : ",Class1," vs ",Class2))+
       theme(plot.title = element_text(size = 16))
    ggsave( image_name , plot1 , width = 10 )
}


use_dat <- subset( dat_diff , class2=="DGC" & class1=="IGC")
image_name <- paste0( out_path , "/IGC_DGC_",subtype,".DiffGene",".pdf" )
diffplot(use_dat = use_dat , image_name = image_name,subtype=subtype,Class2="IGC" , Class1="DGC")


##########################################################################################
##后续直接进行通路富集
result_diff <- dat_diff

computGO <- function(gene_symbol = gene_symbol , type = type){
    gene_df <- bitr(gene_symbol,fromType = 'SYMBOL',toType = 'ENTREZID', 
                    OrgDb = 'org.Hs.eg.db')
    ego <- data.frame(enrichGO(gene_df$ENTREZID,OrgDb = org.Hs.eg.db,ont = 'ALL'))
    ego$type <- type
    return(ego)
}

## 通路富集
enrich_df <- c()
result_diff$CLUSTER <- ifelse(result_diff$log2FoldChange > log2(1.5) & result_diff$padj <0.05 , "DGC_high_expression",
  ifelse(result_diff$log2FoldChange < -log2(1.5) & result_diff$padj <0.05 , "IGC_high_expression","no_sig"))
result_diff <- subset(result_diff,CLUSTER %in% c("DGC_high_expression","IGC_high_expression"))
print(result_diff)
if(nrow(result_diff > 0 )){

      for( clust in unique(result_diff$CLUSTER) ){
        print(clust)
        gene_symbol <- unique(subset( result_diff , CLUSTER == clust )$Hugo_Symbol)
        tmp_go <- computGO(gene_symbol = gene_symbol , type = clust)
        enrich_df <- rbind( enrich_df , tmp_go )
          }

      write.table(enrich_df, paste0(out_path , "enrich.go_filt_",subtype,".txt"), sep = '\t', row.names = FALSE, quote = FALSE)

    ##########################################################################################
    ## 每种类型的通路取前5最显著的
    result <- c()
    top_n <- 5
    for(t_type in unique(enrich_df$type)){
      for( t_path in unique(enrich_df$ONTOLOGY)){
        tmp <- subset( enrich_df , type == t_type & ONTOLOGY == t_path )
        tmp <- tmp[order(tmp$qvalue , decreasing=F),][1:top_n,]
        result <- rbind(result , tmp)
      }
    }

    ## 只看BP通路
    result <- subset(result , ONTOLOGY=="BP")

    gene_ratio <- as.numeric(sapply(strsplit(result$GeneRatio , "/") ,"[" , 1))/as.numeric(sapply(strsplit(result$GeneRatio , "/") ,"[" , 2))
    bg_ratio <- as.numeric(sapply(strsplit(result$BgRatio , "/") ,"[" , 1))/as.numeric(sapply(strsplit(result$BgRatio , "/") ,"[" , 2))
    result$OR <- gene_ratio/bg_ratio
    result$y_p <- -log10(result$p.adjust)
    result <- result[order(result$y_p , decreasing = T),]
    result$Description <- factor( result$Description , levels = unique(result$Description) )
    result$type <- factor(result$type,levels=c("IGC_high_expression","DGC_high_expression"))
    plot <- ggplot(result , mapping = aes(x=Description, y=y_p)) +
        facet_wrap(vars(type) , nrow = 1 , scales = "free_x") +
        geom_bar(stat = 'identity', position = 'dodge') +
        theme_bw()+
        theme(
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.background = element_blank(),
            strip.text = element_text(size = 20),  # 设置分面标签文本大小
            axis.title.x = element_blank(),
            axis.text.x = element_text(color = "black", size = 16, angle = 45, hjust = 1),
            axis.text.y = element_text(color = "black", size = 14),
            axis.title.y = element_text(color = "black", size = 16),
            legend.position = "top",
            legend.title = element_text(color = "black", size = 13),
            legend.text = element_text(color = "black", size = 13))+
              labs(x='GOBP', y='-log10(qvalue)')

    out_name <- paste0(out_path ,  "enrichGO_",subtype,".pdf" )
    ggsave(filename = out_name , plot = plot, width = 15, height = 10)
}  



